TP Simulations de Variables Aléatoires

# pip install plotly

import numpy as np
from matplotlib import pyplot as plt
import numpy.random as nprd
from scipy.stats import norm,cauchy
from IPython.display import clear_output
normal_d,cauchy_d = norm(),cauchy()
from time import time
import seaborn as sns 
import plotly.graph_objects as go
import random as rd

Partie 1

A partir de l'inverse généralisée de la fonction de répartition

L'inverse généralisé de F est donné par : $F^-1 = -\dfrac{1}{\lambda} ln(1-u)$

def simulation_echantillons(n=1000,Lambda=2,norm=1):
    tab_echantillons=[]
    for i in range(n):
        u=np.random.uniform()
        tab_echantillons.append(((-1/Lambda)*np.log(1-u))/norm)
    return tab_echantillons

plt.style.use('ggplot')
particules=simulation_echantillons(1000,2)
sum=np.sum(particules)
particules=[p/sum for p in particules]
plt.hist(particules,bins=100)
plt.title("Echantillons d'une variable aléatoire exponetielle")      

Text(0.5, 1.0, "Echantillons d'une variable aléatoire exponetielle")

norm = lambda a : a/1000
x=np.linspace(0,4,1000)
f=2*np.exp(-2*x)
plt.figure()
plt.plot(x,f,label='densité de probabilité exponentielle pour lambda = 2')
particules=simulation_echantillons(1000,2)
plt.hist(particules,bins=100,density=True)
plt.ylim(0,2)
plt.legend()

<matplotlib.legend.Legend at 0x2d3e7aaba08>

On constate que les échantillons normalisés approchent l'aire sous la courbe de densité. Pour N tend vers l'infini on montre que la somme des aires de l'histogramme tend vers l'aire sous la courbe.
Le problème de cette méthode est qu'elle demande une grande quantité d'achantillons pour estimer l'aire sous la courbe.

Partie 2

M=5
N_sample = 10000
def Accept_Reject(M=5,N_sample = 10000):
    X_gen=[]
    for j in range(N_sample):
        u = nprd.uniform()
        x = nprd.standard_cauchy()
        while 1/M*normal_d.pdf(x)/cauchy_d.pdf(x) < u:
            u=nprd.uniform()
            x=nprd.standard_cauchy()
        X_gen.append(x)
    return(X_gen)



    

X_gen = Accept_Reject()
plt.style.use("ggplot")
plt.hist(X_gen,bins= N_sample//200)
plt.title(" Méthode Accept - Reject")
plt.show()

t=np.linspace(normal_d.ppf(0.0001),normal_d.ppf(0.9999),10000)
X=normal_d.pdf(t)
plt.plot(t,X)
plt.title(" Loi normale centrée réduite")

Text(0.5, 1.0, ' Loi normale centrée réduite')

plt.figure(figsize=[2*5.4, 1*2.8])
plt.hist(X_gen,bins= N_sample//200,color="grey",label="Sample Data",density=True)
plt.plot(t,X,color = "red",label="Densité de $\mathcal{N} (0,1)$",linewidth=2)
plt.title("Estimation de $X \sim \mathcal{N} (0,1)$ à l'aide d'une loi de Cauchy ")
plt.legend()
plt.show()

On constate qu'on approche bien une loi Normale centrée réduite. Ainsi, la méthode $\textit{Accept-reject}$ est une bonne méthode pour générer de nombres aléatoires d'une densité connue à l'aide d'un autre générateur de nombres aléatoires, aussi de densité connue.
On pourra discuter de l'influence de $M$ qui correspond à l'inverse de la probabilité que x soit accepté est égale.

plt.figure()
X_norm=normal_d.rvs(size=N_sample)
plt.hist(X_gen,bins= N_sample//100,color="blue",label="Normal With Cauchy",alpha = 0.4,density=True)
plt.hist(X_norm,bins= N_sample//100,color="red",label="Normal",alpha = 0.3,density=True)
plt.title("Estimation de $X \sim \mathcal{N} (0,1)$ à l'aide d'une loi de Cauchy et \n génération d'une loi normale centrée réduite")
plt.legend()
plt.show()

Qualitativement, on remarque bien que les deux histogrammes se supperposent et on en déduit qu'on a bien un générateur d'une loi normale à partir d'un générateur de nombres aléatoires suivant une loi de Cauchy

Discussion selon le paramètre $M$

m = [1,2,4,8,16]
plt.figure()
T_exec=[]
for _,M in enumerate(m):
    plt.figure(figsize=[2*5.4, 1.5*2.8])
    t_deb=time()
    X_gen = Accept_Reject(M)
    t_exec = time()-t_deb
    plt.hist(X_gen,bins= N_sample//20,color="grey",label="Sample Data",density=True)
    
    T_exec.append(t_exec)
    print(f"Temps d'exécution pour M = {M} vaut {t_exec} secondes" )
    plt.plot(t,X,color = "red",label="Densité de $\mathcal{N} (0,1)$",linewidth=2)
    plt.title(f" Méthode Accept - Reject avec $M={M}$")
    plt.show()
plt.figure()
plt.plot(m,T_exec,color = "black")
plt.xlabel("$M$");plt.ylabel("Temps d'éxécution");plt.title("Temps d'éxécution en fonction de $M$")
plt.semilogx()
plt.show()

Temps d'exécution pour M = 1 vaut 1.6880464553833008 secondes

<Figure size 432x288 with 0 Axes>

Temps d'exécution pour M = 2 vaut 2.4730465412139893 secondes

Temps d'exécution pour M = 4 vaut 5.172055006027222 secondes

Temps d'exécution pour M = 8 vaut 10.372045755386353 secondes

Temps d'exécution pour M = 16 vaut 20.634044885635376 secondes

Plus M est grand, mieux on approxime la loi mais plus le temps d'exécution augmente

3 Algorithme de Box-Muller

Dans l'exercice, au lieu de regarder les résulats à 1 dimension, on utilisera plotly pour regarder les résulats en 2 dimensions.

N_sample = 10000
def Box_Muller(N_sample = 10000):
    X_gen,Y_gen=[],[]
    for j in range(N_sample):
        u1 = nprd.uniform()
        u2=nprd.uniform()
        R= -2*np.log(u1)
        V=2*np.pi*u2
        X_gen.append(np.sqrt(R)*np.cos(V))
        Y_gen.append(np.sqrt(R)*np.sin(V))
    return(X_gen,Y_gen)

X,Y = Box_Muller(10000)

plt.figure()
plt.scatter(X,Y)
plt.show()

plt.figure()
sns.kdeplot(x=X, y=Y,cmap="Blues",fill=True, thresh=0 )
plt.show()

Pour utiliser plotly, installez le à l'aide de $\textit{pip install plotly}$.
Si vous avez besoin de sauvegarder une image, utilisez fig.write_image("Image.png"), mais installez avant kaleido avec $\textit{pip install -U kaleido}$.

h,x1,y1,image = plt.hist2d(X,Y,bins=20,density=True)
fig = go.Figure(data=[go.Surface(z=h, x=x1, y=y1)])
# fig.write_image("Image.png")
fig.show()

plt.imshow(plt.imread("Image.png"))

<matplotlib.image.AxesImage at 0x2d3ecb0ca48>

On remarque qu'on génère bien une gaussienne en 2D à l'aide de la méthode de Box-Muller. Néanmoins, c'est une méthode plus couteuse en temps.

from scipy.stats import multivariate_normal
norm_2d = multivariate_normal([0,0],np.identity(2))
x,y=np.mgrid[-3:3:.1, -3:3:.1]
z= np.dstack((x, y))

fig = go.Figure(data=[go.Surface(z=norm_2d.pdf(z), x=x, y=y,opacity=0.7,colorscale="Greens"),go.Surface(z=h, x=x1, y=y1, opacity=0.4,colorscale='Reds',showscale=False)])
fig.write_image("Image3.png")
fig.show()

plt.imshow(plt.imread("Image2.png"))

<matplotlib.image.AxesImage at 0x2d3ec90a888>

plt.imshow(plt.imread("Image3.png"))

<matplotlib.image.AxesImage at 0x2d3ed0bad48>

La gaussienne 2D se confond avec l'histogramme 2D généré avec la méthode de Box-Muller. En relançant plusieurs fois le programme, on remarque qu'on conserve bien la coincidence entre l'histogramme et la surface de densité de la gaussienne bivariée.

Désormais, générons une loi normal de moyenne -3 et d'écart type 3.

def Box_Muller(N_sample = 10000,mu=-3,sigma=3):
    X_gen,Y_gen=[],[]
    for j in range(N_sample):
        u1 = nprd.uniform()
        u2=nprd.uniform()
        R= -2*np.log(u1)
        V=2*np.pi*u2
        X_gen.append(sigma*np.sqrt(R)*np.cos(V)+mu)
        Y_gen.append(sigma*np.sqrt(R)*np.sin(V)+mu)
    return(X_gen,Y_gen)

X,Y =  Box_Muller(N_sample = 5000,mu=-3,sigma=3)

from scipy.stats import multivariate_normal
norm_2d = multivariate_normal([-3,-3],3**2*np.identity(2))
x,y=np.mgrid[-10:4:.1, -10:4:.1]
z= np.dstack((x, y))
h,x2,y2,image = plt.hist2d(X,Y,bins=30,density=True)

fig = go.Figure(data=[go.Surface(z=norm_2d.pdf(z), x=x, y=y,opacity=0.7,colorscale="Greens"),
                      go.Surface(z=h, x=x2, y=y2, opacity=0.4,colorscale='Reds',showscale=False)])
fig.show()

On obtient ce résulat en remplaçant $\sqrt{R}\sin{V}$ par $\sigma\sqrt{R}\sin{V}+\mu$

4 Algorithme de Generation d'un echantillon normal multivarie de dimension n

from scipy.stats import multivariate_normal

def gen_N_mu_sigma(mu,sigma,d=1000):

    n,X_gen=len(mu),[]

    norm_nd = multivariate_normal(np.zeros(n),np.identity(n))
    A = np.linalg.cholesky(sigma)
    for j in range(d):
        Z = norm_nd.rvs()
        X_gen.append(Z@A+mu)
    return(np.array(X_gen))

import seaborn as sns
X= gen_N_mu_sigma([10,1],[[0.4,0.2],[0.2,0.4]])
plt.figure(figsize=(10,10))
sns.kdeplot(x=X[:,0],y=X[:,1])
plt.show()

mu = np.array([0 ,50 ,100 , 50 ,100 ,200])
sigma = np.array([[11, 10, 5, 9, 4, 2]
        ,[10 ,13, 9, 15, 5, 3]
        ,[5, 9 ,15, 11, 3, 1 ]
        ,[9, 15 ,11 ,21, 6, 4]
        ,[4 ,5 ,3 ,6, 5, 1]
        ,[2 ,3 ,1 ,4, 1, 1]])


X= gen_N_mu_sigma(mu,sigma, 10000)

plt.figure(figsize=(10,5))
sns.kdeplot(X)
plt.show()

for i in range(6):
    print(f"La variance de Z{i} est de {np.diag(sigma)[i]}")
    

La variance de Z0 est de 11
La variance de Z1 est de 13
La variance de Z2 est de 15
La variance de Z3 est de 21
La variance de Z4 est de 5
La variance de Z5 est de 1

Plus la variance est grande, plus la courbe de densité est étalée. C'est ce qu'on remarque sur le graphique ci dessus.

L'amplitude max est donnée par $$\dfrac{1}{\sqrt{2\pi}\sigma}$$
Ainsi, plus sigma est petit, plus l'amplitude max est grande.

import pandas as pd


sns.pairplot(pd.DataFrame(X,columns=[f"Z{i}" for i in range(6)]))
plt.show()

Indice, Corr = [],[]
for i in range(1,len(sigma)):
    for j in range(i):
        Corr.append(sigma[i,j]/(np.sqrt(sigma[i,i])*np.sqrt(sigma[j,j]))) 
        Indice.append(f"{i},{j}")
pd.DataFrame(np.array([np.choose(np.flip(np.argsort(Corr)), np.array(Indice)),np.flip(np.sort(Corr))]).T,columns=["Indices","Corréalation"]).head(5)

	Indices	Corréalation
0	3,1	0.9078412990032038
1	5,3	0.8728715609439696
2	1,0	0.8362420100070909
3	5,1	0.8320502943378437
4	2,1	0.6445033866354896

On peut voir les différentes dispersions selon les diverses sous espaces.
Les variable $Z_3$ et $Z_1$ sont fortement corrélée et on peut le remarquer sur le graphique pairplot ci dessus. On a pu calculer les coefficients de corrélation : $$\rho_{i,j} = \dfrac{\Sigma_{i,j}}{\sigma{i}\sigma{j}}$$

5 Echantillonner suivant une loi de Bernouilli

def Bern(p,d=1000):
    B_gen=[]
    for j in range(d):
        u = nprd.uniform()
        B_gen.append(0 if u <p else 1)
    return(np.array(B_gen))

B = Bern(0.7)
plt.hist(B,density=True)
plt.show()

nb_zeros = (len(B)-np.sum(B))/len(B)
print(f"La fréquence du nombre de zéros est de {nb_zeros}")

La fréquence du nombre de zéros est de 0.718

6 MCMC

t_0 = 500
f = normal_d.pdf
c=0
x_t = [0]
while c <= 20000:
    x,y = x_t[-1],nprd.uniform(-1,1)
    rho = min(f(y)/f(x),1)
    if rd.random() <= rho:
        x_t.append(y)
    else:
        x_t.append(x)
    c+=1
X = x_t[t_0:]
plt.style.use("ggplot")
plt.hist(X,density=True,bins=20)
t=np.linspace(normal_d.ppf(0.0001),normal_d.ppf(0.9999),10000)
X=normal_d.pdf(t)
plt.plot(t,X)
plt.title(" Loi normale centrée réduite")

Text(0.5, 1.0, ' Loi normale centrée réduite')

t_0 = 1000
f = normal_d.pdf
c=0
x_t = [0]
while c <= 20000:
    x,y = x_t[-1],nprd.uniform(-20,20)
    rho = min(f(y)/f(x),1)
    if rd.random() <= rho:
        x_t.append(y)
    else:
        x_t.append(x)
    c+=1
X = x_t[t_0:]
plt.hist(X ,density=True,bins=30)
t=np.linspace(-5,5,10000)
X_f=normal_d.pdf(t)
plt.plot(t,X_f)
plt.title(" Loi normale centrée réduite")

Text(0.5, 1.0, ' Loi normale centrée réduite')

a1 = 1
a2 = 2
mu1 = 10
mu2 = -5
p = 0.3
import scipy.integrate as integrate

f= lambda x: p/(2*a1)*np.exp(-np.abs((x-mu1)/a1)) + (1-p)/(2*a1)*np.exp(-np.abs((x-mu2)/a2))
t_0 = 1000
c=0
x_t = [0]
while c <= 200000:
    x,y = x_t[-1],nprd.normal((mu1+mu2-2)/2,3*(a1**2+a2**2))
    rho = min(f(y)/f(x),1)
    if rd.random() <= rho:
        x_t.append(y)
    else:
        x_t.append(x)
    c+=1
X = x_t[t_0:]
t=np.linspace(-20,20,10000)
X_f=f(t)
plt.hist(X ,bins = 1000,density=True)
a= integrate.quad(f,)
plt.title(" Loi de Laplace")
plt.plot(t,X_f)

[<matplotlib.lines.Line2D at 0x2d3f1749ac8>]

7 Echantillonnage de Monte Carlo parfait

Prenons une variable aléatoire $X$ qui suit une loi normale de moyenne $\mu$ and d'écart type $\sigma$: $$X \sim \mathcal{N}(\mu,\,\sigma^{2})$$ de densité : $p(x) = \frac{1}{\sigma \sqrt {2\pi } }\exp({ - \dfrac{ (x - \mu )^2 } {2\sigma ^2 } })$$

et $$ \mu = \mathbb{E}(X) = \int_\mathbb{R} p(x) dx$$

Ainsi, on a un échantillage de Monte carlo avec $f(x) = x$ .

nb_particule = 100
X = nprd.normal(5,2,size=nb_particule)

mu_hat = np.mean(X)
print(f"La moyenne estimée vaut {mu_hat}")

La moyenne estimée vaut 5.130927466166048

from scipy.stats import norm
nb_particule = 100
V_M = []
for i in range (10000):
    V_M.append( np.mean(nprd.normal(5,2,size=nb_particule)))
plt.hist(V_M,bins=200,density=True)
t=np.linspace(3,7,10000)
X_f = norm(5,2/np.sqrt(nb_particule)).pdf(t)
plt.plot(t,X_f)
plt.show()

La distribution de $V_M$ suit la distribution normal de moyenne $\mu = 5$ et de variance $\dfrac{4}{N}$

nb_particule = 1000
V_M = []
for i in range (10000):
    V_M.append( np.mean(nprd.normal(5,2,size=nb_particule)))
plt.hist(V_M,bins=200,density=True)
t=np.linspace(3,7,10000)
X_f = norm(5,2/np.sqrt(nb_particule)).pdf(t)
plt.plot(t,X_f)
plt.show()

nb_particule = 10000
V_M = []
for i in range (10000):
    V_M.append( np.mean(nprd.normal(5,2,size=nb_particule)))
plt.hist(V_M,bins=200,density=True)
t=np.linspace(3,7,10000)
X_f = norm(5,2/np.sqrt(nb_particule)).pdf(t)
plt.plot(t,X_f)
plt.show()

La variance $\dfrac{4}{N}$ diminue quand le nombre de particule augmente. Ainsi, l'estimation de $\mu$ s'affine quand on prend un plus grand nombre de particule.

8 Echantillonage de Gibbs

import numpy as np
from scipy.stats import norm
from scipy.stats import gamma
import matplotlib.pyplot as plt

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gamma.html

gamma.pdf(2, 2, 0, 1/3)
np.random.seed(123)

m0 = 5
s0 = 4
alpha = 3
beta = 0.5
nbChauffe = 20000
nbEchantillon = 2000 + nbChauffe

m = np.sqrt(s0) * norm.rvs(loc = m0, scale = 1, size = 1)
s = 1 / gamma.rvs(a = alpha, loc = 0, scale = 1 / beta, size = 1)
Y = norm.rvs(loc = m, scale = np.sqrt(s), size = nbEchantillon)

https://numpy.org/doc/stable/reference/generated/numpy.ndarray.flatten.html

plt.hist(Y.flatten(),bins=1000)
plt.grid()
plt.show()

print('m : {}, s = {}'.format(m,s))

m : [8.61458824], s = [0.2877234]

from scipy.io import savemat
mesDonnees1 = {'Y' : Y,
              'm0' : m0,
              's0' : s0,
              'alpha' : alpha,
              'beta' : beta}
savemat('Gibbs.mat', mesDonnees1)
mesDonnees2 = {'m' : m,
               's' : s}
savemat('GibbVraiesValeur.mat', mesDonnees2)

del m0
del s0
del alpha
del beta
del m
del s
del Y

from scipy.io import loadmat
matDonnees = loadmat('Gibbs')
matVraieDonnees = loadmat('GibbVraiesValeur')

Y = matDonnees['Y']
m0 = matDonnees['m0']
s0 = matDonnees['s0']
alpha = matDonnees['alpha']
beta = matDonnees['beta']

m = matVraieDonnees['m']
s = matVraieDonnees['s']
print(len(Y.flatten()))

22000

plt.figure()
plt.hist(Y.flatten(),bins = int(len(Y.flatten())/1000))
plt.grid()
plt.show()

from IPython.display import clear_output

X = np.zeros(shape = (2, nbEchantillon))
X[0, 0] = norm.rvs(loc = m0, scale = np.sqrt(s0), size = 1)
X[1, 0] = 1 / gamma.rvs(a = alpha, loc = 0, scale = 1 / beta, size = 1)
for i in range(1, nbEchantillon, 1):
    M = (s0 / (nbEchantillon * s0 + X[1, i - 1])) * np.sum(Y) + m0 * X[1, i-1] / (nbEchantillon * s0 + X[1,i-1])
    S = s0 * X[1, i-1] / (nbEchantillon * s0 + X[1, i-1])
    X[0, i] = norm.rvs(loc = M, scale = np.sqrt(S), size = 1)
    A = alpha + (nbEchantillon / 2)
    u = (Y - X[0, i])
    v = (Y - X[0, i]).reshape(-1, 1)
    B = beta + 0.5 * (u @ v)
    X[1, i] = 1 / gamma.rvs(a = A, loc = 0, scale = 1 / B, size = 1)

plt.subplot(2,1,1)
plt.plot(X[0, nbChauffe : nbEchantillon])
plt.grid()
plt.subplot(2,1,2)
plt.hist(X[0, nbChauffe : nbEchantillon])
plt.grid()
plt.show()

mHat = np.sum(X[0, nbChauffe : nbEchantillon]) / len(X[0, nbChauffe : nbEchantillon])
print('Estimation de m : ', mHat)

Estimation de m :  8.617519517345205

plt.subplot(2,1,1)
plt.plot(X[1, nbChauffe : nbEchantillon])
plt.grid()
plt.subplot(2,1,2)
plt.hist(X[1, nbChauffe : nbEchantillon])
plt.grid()
plt.show()

sHat = np.sum(X[1, nbChauffe : nbEchantillon]) / len(X[1, nbChauffe : nbEchantillon])
print('Estimation de s : ', sHat)

Estimation de s :  0.28727648846856063